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Abstract 



Wave packet propagation in the basis of interpolating scaling functions (ISF) is stud- 
ied. The ISF are well known in the multiresolution analysis based on spline biorthog- 
onal wavelets. The ISF form a cardinal basis set corresponding to an equidistantly 
spaced grid. They have compact support of the size determined by the order of the 
underlying interpolating polynomial that is used to generate ISF. In this basis the po- 
tential energy matrix is diagonal and the kinetic energy matrix is sparse and, in the ID 
case, has a band-diagonal structure. An important feature of the basis is that matrix 
elements of a Hamiltonian are exactly computed by means of simple algebraic trans- 
formations efficiently implemented numerically. Therefore the number of grid points 
and the order of the underlying interpolating polynomial can easily be varied allowing 
one to approach the accuracy of pseudospectral methods in a regular manner, simi- 
lar to high order finite difference methods. The results of numerical simulation of a 
H+H2 collinear collision show that the ISF provide one with an accurate and efficient 
representation for use in the wave packet propagation method. 
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1 Introduction 



A direct on-the-grid solution of the time-dependent Schrodinger equation (TDSE) has be- 
come a common tool in quantum chemistry. In a dynamical context, it provides with quan- 
titative predictions on the efficiency of different pathways of chemical reactions and deepens 
our understanding of their details [U El El- Approaches based on the TDSE in connec- 
tion with the filter diagonalization technique jlj are also used in the static context in order 
to compute states of complex molecules Ej. To simulate the time evolution of a wave 
packet, repeated computations of the action of the system Hamiltonian H, or its exponen- 
tial exp(— iAtH), on the wave function ^ are required. Here At is a time step. Therefore a 
lot of effort has been devoted to developing accurate and efficient methods to reduce com- 
putational costs of computing H^fr. The numerical techniques used can be classified into 
several categories: (i) finite differences (FD), (ii) finite elements (FE), (iii) pseudospectral 
global grid representation approaches such as the discrete variable representation (DVR) 
and Fourier grid Hamiltonian method (FGH) El El E3 E] • Finite differences and finite 
elements lead to a sparse Hamiltonian matrix but exhibit a slow algebraic convergence with 
the number of grid points. The pseudospectral approaches result in dense Hamiltonian ma- 
trices, which, a priori, increases the number of operations needed to compute ff\P. At the 
same time, the exponential convergence with the number of grid points counterbalances the 
aforementioned drawback. This is why pseudospectral methods are so widely used in time- 
dependent molecular dynamics [H El U2j, as well as in stationary S-matrix [T31 EH HHj 
and eigenvalue |6 calculations. In particular, the FGH method based on the fast Fourier 
transform (FFT) algorithm is very advantageous since for a mesh of iV points the action of 
the kinetic energy operator is computed by N log 2 N elementary multiplications and is easily 
implemented numerically [Tfi] . 

When discussing the slow algebraic convergence of finite differences and finite elements, 
one usually refers to the convergence with the number of grid points while the order of 
finite differences or the order of the polynomial for finite elements is kept fixed. In fact, an 
exponential (spectral) convergence can be achieved if not only the number of grid points, but 
also the order of the underlying polynomial is increased 0[T7j. Exploiting this property, as 
well as the sparsity of the Hamiltonian matrix and the possibility to distribute mesh points 
only along the reaction path, thus reducing the size of calculation, makes the high order FD 
and FE a competitive alternative to pseudospectral methods. Much work has recently been 
devoted to developing these techniques [TBI 113 1201 [2H 123 123 12H 

In this paper we present a treatment of reactive scattering based on the wave function 
representation in the basis of interpolating scaling functions corresponding to interpolating 
(spline) wavelets j2H] • Basis functions are generated from a single function, called the scaling 
function, by appropriate scalings and shifts of its argument. The scaling function is a solution 
of a functional equation which is found iteratively by using interpolating polynomials of a 
specific order. By construction, the basis functions have compact support with the width 
determined by the order of the interpolating polynomial and the resolution level (the mesh 



2 



step). The basis is biorthogonal and cardinal, which leads to a simple representation of the 
wave function. The resulting Hamiltonian matrix is sparse. The potential energy matrix 
appears to be diagonal, while, in the one dimensional case, the kinetic energy matrix is band 
diagonal with the band width determined by the size of support of the scaling function. 
Matrix elements of the kinetic energy operator can be evaluated exactly via simple algebraic 
operations so that the order of the underlying interpolating polynomial and the step of the 
mesh can easily be varied allowing one to achieve a fast convergence. Thus our algorithm 
offers flexibility similar to that achieved in the case of finite differences of an arbitrary 
order with the Fornberg algorithm |2H]. In fact, the proposed approach can be seen as 
an alternative to high-order finite difference techniques. Applications of the method are 
illustrated with the example of a collinear H+H2 collision. 



2 Theory 

2.1 Biorthogonal spline bases 

Here we summarize some necessary facts about interpolating biorthogonal bases of scaling 
functions in the space of square integrable functions and describe the algorithm to construct 
such bases. Biorthonormal bases of scaling functions are used in the multiresolution analysis 
associated with interpolating wavelets. A more detailed description of such bases can be 
found in the mathematical literature, e.g., [23 123 EH] • 

In general, a biorthogonal basis consists of two sets of elements <j) a (x) and (j) a (x) where 
the index a labels the basis elements. Any function can be decomposed into a linear 
combination of the basis functions 4> a {x), 

=Y. S a<l>a(x) , (1) 
a 

where the decomposition coefficients are determined by the dual basis functions 

s a = J dxMx)^(x) . (2) 
The basis is (bi)orthogonal in the sense that 

dx 4> a (x)(f) b (x) = 5 ab . (3) 



Consider a special class of biorthogonal bases that provide a multiresolution analysis. A 
biorthogonal basis with a multiresolution analysis is generated by a scaling function <f>{x) 
and its dual <j>(x) which satisfy accordingly the equations (the scaling relations) 

<t>(x) = 2j2h k ^(2x-k) , (4) 

k 

4>{x) = 2£M(2x-A;) , (5) 

k 
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where the real coefficients hk {hk) are called a (dual) filter. The scaling function and its dual 
are required to be orthogonal in the sense that 

dx 4>(x)(f>(x - j) = 5 0j (6) 

for all integers j. The orthogonality relation imposes a condition on the filters which is 
readily deduced from (0) by substituting the scaling relations (J1J) and (jSJ) are rescaling the 
integraion variable. 

Consider two sets of functions, labeled by two integers n and j, 

<j> nd {x) = T'^iTx - 3 ) , 4> nj (x) = 2 n / 2 4>(2 n x - j) . (7) 

Subspaces V n of the space of square integrable functions spanned by 0„j with a fixed value 
of n form a ladder structure 

• • • C V n C V n+1 C • • • . (8) 
It is straightforward to convince oneself that <p n j form a biorthogonal basis in V n : 

dx(f) n j(x)(f) n j / (x) = Sjji . (9) 

Any function \^ can be projected on a subspace V n , 

P n : $?(x) ^n(x) = ^2 s n:j (f) n j(x) , s n j = / dx <j> n j (x)^f(x) . (10) 

3 J 

Taking successively larger values of n allows one to reproduce a successively finer structure 
of ty. Thus, the index n specifies a resolution level. If the filters are finite, then the scaling 
function has finite support and so do the basis functions <p n j. The index j is then naturally 
associated with the position of support of <p n j. 

A special class of biorthogonal bases is obtained when filters are finite and of a special 
form 

hk = ^(fc/2), fc = 0,±l,±2,...,±m J (11) 
h k = 4o , <j>(x) =S(x) . (12) 
The expansion coefficients are simply values of the function at dyadic lattice sites 



s 



2 -n/2^( 2 -™7) . (13) 



Larger values of the resolution level n correspond to finer grids. In what follows m is chosen 
to be odd for convenience. 

To find the filter hk and an explicit form of the scaling function 0(x), Eq. (jlj) is solved 
iteratively for <p{x). First, one observes that the equation is satisfied at the integer valued 
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argument by (f)(1) = 5 i with ho = ~. Thus, the scaling function vanishes at integral x except 
for x = where 0(0) = 1 as is required by the orthogonality condition © with our choice of 
the filters (fTTj) and (fT2j) . To compute (p(x) at half-integer values of the argument, 4>(j + 1/2) 
(j is a fixed integer), one uses the polynomial (spline) interpolation with polynomials P£j/ of 
order m passing through M = [m + 1)/2 integer points neighboring to x = j + 1/2 to the left 
and M integer points to the right: x = j — M + l,j — M + 2, j + M. Having found such 
a polynomial, we set <p(j + 1/2) = P^)(j + 1/2). Note that for different points x = j + 1/2, 
polynomials P$ are different. To find 0(j±l/4), the same procedure is applied, but now the 
values of at M half-integer points neighboring to x — j ±1/4 to the left and M half- integer 
points to the right are used to construct the corresponding interpolating polynomial. For 
example, for x — j + 1/4 the sequence will be: x = j — (M — l)/2, j — (M — 2)/2, j + 
M/2. In other words, any x can be squeezed into successively smaller intervals of length 
1/2, 1/4, 1/2^, N — > oo. The limiting procedure allows one to compute <f>(x) at any x, 
in principle. 

An important property of <p(x) is that it has compact support, an interval of width 
D = 2m, that is, 4>(x) = for all \x\ > m as can be inferred from the construction procedure 
described above. 

In Fig. 1 we show an example of the scaling function 4>(x) = 0o,o( a; ) where 4>o,o{ x ) is 
the basis function centered at position x = and corresponding to the zero resolution level 
(the mesh mesh equals 1). The length of the filter is m = 15 so that support of the scaling 
function is the interval —15 < x < 15. 



2.2 Hamiltonian matrix 

Consider first a simple ID example. Let if be a Hamiltonian of a system in the coordinate 
representation. A solution of the Schrodinger equation id/dt *f>(x,t) = H^f(x,t) for a given 
initial wave packet ty(x,t = 0) is approximated by its projection into a finite dimensional 
subspace spanned by n j where j enumerates basis functions whose support lies in a box, x G 
[0, L], that is, *&(x, t) = J2j s n ,j(t)4>n,t( x )- Each basis function <p n j has support of the length 
D„ = 2m/2 n . Therefore the number of the basis functions N is given by: iV = 2 n L/2m. 
Our choice of the basis implies also zero boundary conditions for the wave function. The 
initial wave packet *&(x, t = 0) is projected into the corresponding subspace of V n according 
to fTU|) to determine s„j(0). The Hamiltonian is projected by the rule H P n HP n and 
becomes a finite matrix with elements 

H$ = Jdx^ n>k (x)H<f) nj (x) . (14) 

The Hamiltonian matrix acts in a vector space of the expansion coefficient s n j(t). Solving 
the time dependent Schrodinger equation thus implies computing s n (t) = exp 
where s n is regarded as a vector with components s n j and as a matrix with elements 
given in (JTJJ). 
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A typical Hamiltonian is a sum of kinetic and potential energies. For matrix elements of 
the potential energy V we have 

V k f = J dx^{x)V{x)4> n>j {x) = V(2~ n j)S kj . (15) 

This is a great advantage of the basis under consideration: The potential energy matrix is 
diagonal. To compute the kinetic energy matrix we apply a general procedure to compute 
derivative operators d l /dx l in a basis of compactly supported (spline) wavelets jH0j- By 
rescaling the integration variable, and using Eqs. (JZJ) and l|12j) one can find 

D$f = [ dxj) njk (x)(d/dx) l (f) nd (x) = 2 nl J dx4>{x){d/dx) l (f){x-k + j) 

= 2 n \d/dx ) l <j>(x -k + j) \ x=0 = 2 nl D% . (16) 
Using the scaling relation for the scaling function (jlj) we infer 

Df = ^j^ej:^, (I?) 

k j 

4 } = 2 l ^h^ = 2 l <j>{i - j/2) . (18) 

Thus, DW is an eigenvector of the matrix corresponding to the eigenvalue 1. 

Finally, we need a normalization of the vector D^ l \ Note that, since the value of <p a t 
any x is given by a polynomial of order m, the monomial x l , I < m, should be a linear 
combination of <fio,j( x ) — 4>{ x ~ j); that is, x l = J2j s oj4>o,j( x ) with soj = j 1 . Differentiating 
this relation / times, multiplying by <p(x) and integrating over x we obtain the normalization 
relation 

l\ = J23 l Df ] . (19) 

j 

The matrix D^ n ' 1 ^ satisfies the symmetry relation = (— l) l D^ l \ which follows from 

(116(1 after changing the integration variable x — > — x and making use of <p(x) = <p{—x) (the 
same for the dual). In particular, D^' 2 ^ is a symmetric matrix. Since support of 4>{x) lies 
within \x\ < m, the D^' 1 ^ matrix is band diagonal: D^p = 0, if \k — j\ > m as follows from 
Eq. ((Tfij) . Therefore the action of the Hamiltonian of the system = —^D^ n ' 2 ^ + on 
the wave function s n requires N x 2m elementary multiplications. The above approach allows 
easy and fast evaluation of the Hamiltonian matrix for any desirable resolution level n (any 
number of the basis functions) and any length of the filter m (any interpolating polynomial 
order) . 

A multidimensional generalization is obtained by taking the direct product of raj for 
every independent variable Xj. The resolution level n« may be chosen independently for 
every variable For example, in the two dimensional case, the basis consists of func- 
tions 4>n 1 ,j 1 (xi)(f) n2t j 2 (x2). The spline order m may also be chosen independently for every 
coordinate. In other words, the scaling functions for each Xi may be different. 
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3 Numerical results 



3.1 Harmonic oscillator example 

In Fig. 2 we show the results of a test calculation of 30 first eigenvalues of a ID harmonic 
oscillator, H = — + \- We use a mesh in the interval — L < x < L with L = 10 a where 
cio is the Bohr radius. The Hamiltonian matrix in the basis of interpolating scaling functions 
4> n ,j(x) has been calculated according to the procedure described in the previous section 
and diagonalized yielding the eigenstate energies. The results of the present approach are 
compared with those obtained by the Fourier Grid Hamiltonian (FGH) approach. In the FGH 
method, the convergence is reached with 80 points of the grid and the error of the eigenvalue 
calculation is basically determined by the precision of the diagonalization procedure. In the 
ISF basis and low order m, convergence with the number of basis functions iV (N = 2 n x 2L, 
where n is the resolution level) is slow. At the same time the convergence can dramatically 
be improved by increasing the order m of the interpolating polynomial, i.e. by increasing the 
band width of the band-diagonal kinetic energy matrix. This observation is in line with the 
results reported by several authors for finite differences and finite elements where it has been 
shown that the pseudospectral convergence can be reached by increasing the order of finite 
differences or the order of an underlying polynomial in finite elements [°Q [""7J EH • 
Thanks to the biorthogonality of the basis in our case, increasing the order m and the 
resolution level n is a quite simple procedure given by Eqs. (JT3J) (JTHJ) which is easy to 
implement numerically. 

In the case of the harmonic oscillator the FGH method outperforms the ISF method. The 
latter, as well as finite difference methods, might compete with pseudospectral approaches 
for calculations involving complex reaction paths. Indeed, the FGH method, for example, 
generally uses hypercubic grid domains. In the present approach the calculation volume 
can sufficiently be reduced by the choice of the basis functions 4>n,j whose support lies in the 
vicinity of the reaction path, which is again a simple procedure because of the biorthogonality 
of the basis. We now turn to one such example. 



3.2 H+H2(f = 0) collinear collision 

Here we present a wave packet propagation treatment of an H+H 2 collision in collinear 
geometry with the energy of the collisional Hydrogen atom between 0.2 eV and 1.1 eV. The 
system is described by Jacobi coordinates with r± being the distance between the two 
Hydrogen atoms in the molecule and r 2 being the distance of the collisional Hydrogen atom 
from the molecular center of mass. We use the Lanczos method [3T1 13*2*] to solve the time 
dependent Schrodinger equation (atomic units are used) 

+ V(r u r 2 ))*(r ll r 2 ;t) , (20) 



2m 



H 



39^ 
+ 2dr\ 
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where m# is the mass of the Hydrogen atom. The interaction potential V A (r 1 ,r 2 ) is taken 
from [33]. A typical size of the mesh is 32 x 32 a . An absorbing potential |Ml 133] is 
introduced at the grid boundaries for large r± and r 2 to avoid interference of the simulated 
wave packet with its reflection from the grid boundary. The initial state corresponds to the 
H2 molecule in the v = vibrational state and an impinging Gaussian wave packet in the 
reaction channel (r 2 ). 

In the FGH method, convergent results are obtained with 192 points in the T\ coordinate 
and 256 points in the r 2 coordinate. Two types of calculation have been performed by 
the ISF method. The first one corresponds to a rectangular grid where the basis functions 
are chosen as $„j lj2 = 0n J1 (ri)0 ni j 2 (r 2 ). The resolution level n = 4 has been used, which 
corresponds to the mesh step 1/2 4 = 0.0625 a so that in total there are 360 x 512 basis 
functions. The order of the underlying interpolating polynomial (size of the filter) has been 
set to m = 21. In the second simulation the ISF basis has been chosen in such a way that 
the basis functions have support in the potential energy valley inside the 4.6 eV potential 
energy level curve as shown in Fig. 3. In this case only 48600 basis functions are needed. 
As a result, the computational time has been reduced by 3.8 times, which brings it to the 
level comparable to the FGH treatment. Results are presented in Fig. 4. We find that with 
our choice of the basis both simulations based on the interpolating scaling functions yield 
the total reaction probability R which nicely agrees with the FGH method. The absolute 



difference 



Rfgh—Risf 
Rfgh 



is better than 0.1%. 



4 Conclusions 



We have shown that the wave function representation in bases of interpolating scaling func- 
tions can efficiently be used in wave packet propagation studies of molecular dynamics. The 
ISF originate from a multiresolution analysis associated with interpolating (spline) wavelets 
[2*7| |2Hl EH! The ISF have a finite support whose width is determined by the order of the 
underlying interpolating polynomial and the resolution level (the mesh step). The basis is 
cardinal, which leads to a simple representation of the wave function. The resulting Hamil- 
tonian matrix is sparse. In particular, the potential energy matrix is diagonal, while for 
the kinetic energy matrix there is an efficient and simple algebraic procedure for its exact 
computation so that the order of the underlying interpolating polynomial and the mesh step 
can easily be varied to achieve fast convergence. Thus, our algorithm provides one with 
flexibility similar to that achieved in the case of finite differences of an arbitrary order by, 
e.g., the Fornberg algorithm |2E]- In the ID case the kinetic energy matrix has a band 
diagonal structure with the width of the band of off-diagonal elements given by the size of 
support of the scaling function. In our example of a collinear H+H 2 collision, a considerable 
reduction of computational costs has been reached when the basis functions are chosen in 
such a way that their support is localized only in the vicinity of the reaction path. The very 
possibility of the latter procedure and the simplicity of its numerical implementation, thanks 
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to the biorthogonality of the basis, is a general feature of the method proposed which can be 
used in simulations of systems with more complex reaction paths to reduce computational 
costs. 
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Figure captions 

Fig. 1. Scaling function (f>(x) for the filter length m = 15. 

Fig. 2 A ID harmonic oscillator problem. A relative error in the eigenvalue calcula- 
tion as a function of the oscillator quantum number. Lines: The FGH method with different 
numbers A" of mesh points; the gray line is for A" = 40, the dashed line for A" = 80, and 
the black line for A" = 160. Symbols and lines with symbols: The ISF method. Circles 
show the results obtained with n = 2 resolution level corresponding to the A^ = 80 scaling 
functions basis (N = 80 points of the mesh). Triangles show the results obtained with n = 3 
resolution level corresponding to the A" = 160 scaling functions basis (A" = 160 points of the 
mesh). Gray solid symbols: Results obtained with the filter length m = 7 (the interpolating 
polynomial order); open symbols: Results obtained with the filter length m — 15; black solid 
symbols: Results obtained with the filter length m — 21. The line with black solid symbols: 
results obtained with the filter length m — 41. 

Fig. 3 A schematic representation of the arrangement of of the mesh (positions of 
the basis functions) in the wave packet propagation treatment of a collinear H-H 2 collision. 

Fig. 4. Reaction probability calculated by different methods for the H+H 2 (t> = 0) 
collinear collision as a function of energy. Solid line: Fourier Grid Hamiltonian method; 
Open circles: Results obtained by ISF method in the rectangular mesh; Triangles: Results 
obtained by the ISF method in the mesh arranged along the reaction coordinate (path) as 
depicted in Fig. 3. 
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